% This script rescales and plots impulse response functions for 
% contemporaneous or news TFP shocks.

% Christoph Gortz


close all
clear all

%% Initialisation

 
shockPlot = 'z';             % Shock process to be plotted
exo_shockPlot = 'ez8';        % Exogenous shock to be plotted
shockOther = 'v';             % Other Shock process with trend
exo_shockOther = 'ev';        % Other Exogenous shock with trend
vars = [                    % Variables to be plotted
    'y         ';
    'c         ';
    'r         ';
    'pi        ';
%     'sc     ';
%     'si     ';
    'spreadobsc';
    'spreadobsi';
    %
    'inve      ';
    'invec     ';
    'invei     ';
%     'nc     ';
%     'ni     ';
    'qc        ';
    %
    'pinfc     ';
    'pinfi     ';
    'lab       ';
    'labc      ';
    'labi      ';
%     'lrc    ';
%     'lri    ';
%     'spkc   ';
    %
%     'w      ';
    'qi        ';
%     'rrate     ';
    'mcc       ';
    'mci       ';
    ];
nVars = size(vars,1);
IRFtitle = [                % Title of figures corresponding to variables
    'Output                   ';
    'Consumption              ';
    'Nom. Interest Rate       ';
%     'C-Sector Financial Claims';
%     'I-Sector Financial Claims';
    'Rel. Price of Investment ';
    'C-Sector Spread          ';
    'I-Sector Spread          ';
    %
    'Total Investment         ';
    'C-Sector Investment      ';
    'I-Sector Investment      ';
%     'C-Sector Bank Capital    ';
%     'I-Sector Bank Capital    ';
    'C-Sector Price of Capital';
    %
    'C-Sector Inflation       ';
    'I-Sector Inflation       ';
    'Total Hours              ';
    'C-Sector Hours           ';
    'I-Sector Hours           ';
%     'C-Sector Leverage Ratio  ';
%     'I-Sector Leverage Ratio  ';
%     'C-Sector Valuation Shock ';
    %
%     'Real Wage                ';
%     'C-Sector Price of Capital';
    'I-Sector Price of Capital';
%     'Real Interest Rate       ';
    'C-Sector Price Markup    ';
    'I-Sector Price Markup    ';
    ];



% Names of .mat files with IRFs to be plotted in same figure 
dataName = ['Output_Baseline3              '  
            'Output_noFIestimated_sameStd_2']; % estimated model with no financial sector -- same std as baseline
figName = ['IRF_baseline_vs_noFIestimated3_v2' , exo_shockPlot];



figDimensions = [3, 6];     % Dimensions of figure - needs to be consistent
                            % with nVars
irfLength = 40;             % Length of generated IRFs 


%% LOAD INFO AND DETREND
for iData = 1:size(dataName,1)
    eval([ 'load  ' dataName(iData,:) ';' ])
    
    

% Load variables' irf data into cell
for iVars = 1:nVars
    if iData == 1
    eval([ 'irfMean(iVars,1) = {oo_.irfs.' strcat(vars(iVars,:)) '_' exo_shockPlot '};' ])
    elseif iData == 2
        if iVars == 5 
            irfMean(iVars,1) = { zeros(1, 40) };
        elseif iVars  == 6
            irfMean(iVars,1) = { zeros(1, 40) };
        else
            eval([ 'irfMean(iVars,1) = {oo_.irfs.' strcat(vars(iVars,:)) '_' exo_shockPlot '};' ])
        end
    end
end

for i = 1:numel(irfMean);
    Variables{i,1} = [irfMean{i,1}];
end
% Load shock's irf data
eval([ shockPlot 'Mean = oo_.irfs.' shockPlot '_' exo_shockPlot ';' ])
eval([ shockOther 'Mean = oo_.irfs.' shockOther '_' exo_shockPlot ';' ])



% Transform series of shock z to isolate zd:
vd =[]; vp = 0;
zd =[]; zp = 0;

for m = 1:numel(zMean)
    % Change of v and z in percent
    vp = vMean(1,m)+vp;
    vd = [vd; vp];
    zp = zMean(1,m)+zp;
    zd = [zd; zp];
end    
vd = vd';
zd = zd';


% Descale stationarised variables: K_t = K_t^~*Z_t
yd = Variables{1,1} + (alfac/(1-alfai))*vd + zd;
cd = Variables{2,1} + (alfac/(1-alfai))*vd + zd;
rd = Variables{3,1};
pid = Variables{4,1} + zd - ((1-alfac)/(1-alfai))*vd;
spreadcd = Variables{5,1};
spreadid = Variables{6,1};
inved = Variables{7,1} + (1/(1-alfai))*vd;
invecd = Variables{8,1} + (1/(1-alfai))*vd;
inveid = Variables{9,1} + (1/(1-alfai))*vd;
qcd = Variables{10,1} - ((1-alfac)/(1-alfai))*vd + zd;
pinfcd = Variables{11,1};
pinfid = Variables{12,1};
labd = Variables{13,1};
labcd = Variables{14,1};
labid = Variables{15,1};
qid = Variables{16,1} - ((1-alfac)/(1-alfai))*vd + zd;
% rrated = Variables{17,1};
% zcapcd = Variables{12,1};
% zcapid = Variables{13,1};
% rkcd = Variables{14,1} + zd - ((1-alfac)/(1-alfai))*vd; % Is this correct?
% rkid = Variables{15,1} + zd - ((1-alfac)/(1-alfai))*vd;
mccd = Variables{17,1} * (-1);                                 % Is this correct?
mcid = Variables{18,1} * (-1);
% wd = Variables{18,1} + zd + alfac/(1-alfai)*vd;
% lrcd = Variables{21,1};
% lrid = Variables{22,1};
% ncd = Variables{23,1} + (alfac/(1-alfai))*vd + zd;
% nid = Variables{24,1} + (alfac/(1-alfai))*vd + zd;


% Final Detrended Time Series of IRFs
eval([ 'VariablesDetrended_' dataName(iData,:) '= [yd; cd; rd; pid; spreadcd; spreadid; inved; invecd; inveid; qcd; pinfcd; pinfid; labd; labcd; labid; qid; mccd; mcid;];' ]);

end
%% Plot Detrended IRFs

%% Plot IRFs
figure1 = figure('Color',[1 1 1]); % Generates white background around plot

    for iVars = 1:nVars
        varStr = strcat(vars(iVars,:));
        plotVar = strcat(vars(iVars,:));
        irfTitle = strcat(IRFtitle(iVars,:));
        % Subplots
        subplot(figDimensions(1),figDimensions(2),iVars);
        for iData = 1:size(dataName,1) % Loop over different datasets of model versions
        eval([ 'irfData = VariablesDetrended_' dataName(1,:) '(iVars,:);' ]);
        eval([ 'plot(irfData,''-k'',''linewidth'',1.5);' ]);
        hold on
        eval([ 'irfData = VariablesDetrended_' dataName(2,:) '(iVars,:);' ]);
        eval([ 'plot(irfData,''-or'',''linewidth'',1.5);' ]);
        if iData == 3
            eval([ 'irfData = VariablesDetrended_' dataName(3,:) '(iVars,:);' ]);
            eval([ 'plot(irfData,''--b'',''linewidth'',1.5);' ]);
        end
        end
        axis tight
        hold off
        eval( 'title( irfTitle );' );
    end
set(gcf, 'Position', get(0,'Screensize')); % Maximise Figure


